Stable and unstable vortices in multicomponent Bose-Einstein condensates 
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We study the stability and dynamics of vortices in two-species condensates as prepared in the recent 
JILA experiment ( M. R. Andrews et al, Phys. Rev. Lett. 83 (1999) 2498). We find that of the 
two available configurations, in which one specie has vorticity m = 1 and the other one has m = 0, 
only one is linearly stable, which agrees with the experimental results. However, it is found that 
in the unstable case the vortex is not destroyed by the instability, but may be transfered from one 
specie to the other or display complex spatiotemporal dynamics. 
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Vortices appear in many different physical contexts rang- 
ing from classical phenomena such as fluid mechanics |lj and 
nonlinear Optics M to purely quantum phenomena such as 
superconductivity and superfluidity [g]. In the last two years 
more than 100 papers concerning vortices have been published 
in Physical Review Letters, which is a naive way to appreciate 
the importance of this subject in Physics. 

A vortex is the simplest topological defect one can con- 
struct: in a closed path around a vortex, the phase of the 
field undergoes a 27r winding and stabilizes a zero value of the 
field placed in the vortex core. The vortex is stable because 
of topological constraints; removing the phase singularity im- 
plies an effect on the boundaries of the system which cannot 
be achieved using local perturbations. 

Vortices are central to our understanding of superfiuidity 
and quantized fiow. This is why after the experimental re- 
alization of Bose-Einstein condensates (BEG) with ultracold 
atomic gases n the question of whether atomic BEG's are 
superfluids has triggered the analysis of vortices. The main 
goals up to now have been to propose a robust mechanism 
to generate pH and detect vortices Q. But another impor- 
tant research area is the analysis of vortex stability |8|-[10[, to 
which this works contributes. 

Although most of the theoretical effort concerning vortices 
has been focused on single component condensates, the first 
experimental production of vortices in a BEG was attained 
using a two-species *^Rb condensate |ll|]. In this experiment 
the condensed cloud is made up of atoms in two different 
hyperfine levels, denoted by |1) and |2). Since the scattering 
lengths are different both states are not equivalent. As a 
consequence, while each specie can host a vortex, in Ref. Illl] 
it was shown that a single vortex is stable only when it is 
placed in the component with the largest scattering length, 
|1). The other possibility, which has the vortex in |2), leads 
to some kind of instability. 

Our intention is to prove in this paper that the instabil- 
ity is purely dynamical and can be explained with a mean 
field model which does not include any type or dissipation. 
We will achieve this goal in three major steps. First, we pro- 
pose a model based on coupled Gross-Pitaevskii equations and 
solve the stationary equations in two and three-dimensional 
setups. We obtain the lowest energy stationary states that can 
be qualitatively identified with the ground state and the two 
realizations of Ref. |ll| . Next we study the stability of each 
state under small perturbations using linear perturbation the- 



ory. Our main result is that only the experimentally stable 
configuration is also linearly stable. Finally, using numeri- 
cal simulations, we study more realistic conditions in which 
the condensate suffers moderate to strong perturbations. We 
show that there is a good agreement with experiment and also 
that the dynamics is very rich and depends on the dimension- 
ality of the system and the intensity of the perturbations. 

The model- In this work we will use the zero temperature 
approximation, in which collisions between the condensed and 
non condensed atomic clouds are neglected. In the two species 
case this leads to a pair of coupled Gross-Pitaevskii equations 
(GPE) 
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are the corresponding scatter- 
ing length. To simplify the formalism and in analogy with 
experiments we assume that both potentials are spherically 
symmetric, i. e. Vi{r) — V2(r) = ^muj'^{r^ + z^). 

Following the experiments, we will present results for equal 
number of particle on each specie, A'^i = N2 = A*', which 
translates to 



\M^fdr= / IM^l^'df 



(2) 



after a proper rescaling of tpi and ip2 by A^. To simplify 
the analysis we change to a new set of units based on 
the trap characteristic length, ao — yjh/rnuj, and period, 
T — 1/uj. In this set of units the nonlinear coefficients are 
Uij = AixOij ^J NiNj /ao . For the JILA experiment, in which 
u) = 2n X 7.8 ± O.lHz, we have that r ~ 20.4 ms is the new 
unit of time. 

The whole study, including the linear stability analysis and 
the numerical simulations, was performed in two- and three- 
dimensional systems. We have studied the system up to values 
of Uij < 5000, which are of the order of magnitude of the ex- 
periments. Nevertheless, the linear stability analysis and the 
simulations change little for the strongest interactions and we 
expect our results to be still valid for a larger number of par- 
ticles {N ^ 10^). Regarding the intensity of the nonlinearity, 
we have used the scattering lengths of *^Rb which appear in 



Ref. [[L2| . These values give us a precise line in the parameter 
space 



U^g 
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It is remarkable that because of the relation un > ui2 > U22 
the experiment is performed in a regime in which the first 
component chases the second one, which rejects mixing with 
the chaser. This means that a "desired" configuration has the 
first component spread over the largest part of the space. As 
we will see, this has important consequences for the dynamic 
of the states. 

Search of solutions.- We are interested in stationary con- 
figurations in which each component has a well defined value 
of the angular momentum. Such states have the time and 
angular dependence factored out 
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These functions satisfy a nonlinear set of coupled PDE 
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with i — 1,2. Our focus will be on three particular con- 
figurations, which are the lowest energy states with vortici- 
ties (771,1,7112) = (0, 0), (1,0), (0, 1). They correspond to the 
ground state of the double condensate, and to the single vor- 
tex states for the |1) and |2) species, respectively. 

To find the radial and longitudinal dependences of the 
wave functions, (f>i{r,z), we expand them approximately on 
a finite subset of the harmonic oscillator basis, (f)i{r,z) 
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states for each of the (7111,7112) pairs of vorticities. The de- 
tails of the method applied to single vortex systems can be 
found in Ref. [|10|. As a result one obtains the desired eigen- 
functions, chemical potential and energies [Fig. Wa)] of each 
(7iil,7n2) ground state. 

Linear stability analysis.- Let us study the behavior of the 
stationary states under infinitesimal perturbations (e.g. in 
the initial data, small amounts of noise, etc.). To do so, we 
define the excitations as 



Mr, z, 6) = 6-'"'*+""'" (<^,(r, z) + e-'^*+'"«a.(r, z)) , (6a) 
Mr, z, 9) = e'"'*"""'» (0,(r, z) + e~'^'''"'Mr, z)) . (6b) 



Then we introduce Eqs. (m) into Eq. (1) and keep the 
0{a) and C(/3) terms. In the end we reach an equation 
for the perturbations W — [a, a*,/3, /3*] which is of the type 
iJdtW = HW. Here H is an hermitian operator and ij is an 
anti-hermitian operator. We will search a Jordan basis such 
that XW = —iJHW. The lack of such a Jordan form leads to 
polynomial instabilities, while the presence of non-real eigen- 
values is a signal of the exponentially growing instabilities. 
All other modes give us the frequencies of the linear response 
of the system to small perturbations. 

This analysis is formally equivalent to the Bogoliubov sta- 
bility analysis of the states under consideration. However, as 
it was shown in [ho|, our perturbation is of order 0{a^,P^) 
in the energy and thus the exponentially unstable modes lay 



along directions of constant energy and thus a Bogoliubov ex- 
pansion of the Hamiltonian cannot account for the instability, 
which can only be reached by working directly on the GPE. It 
is thus important to confirm our linear stability results with 
numerical simulations of the whole system. 

Linear stability results.- When we perform the diagonaliza- 
tion on (7111 , 7712) = (0,0) we obtain that all of the eigenvalues 
are positive numbers, which implies that (0, 0) is at least a 
local minimum of our energy functional -furthermore, it is 
the ground state of the system and thus globally stable. 

Next we studied the (1,0) family [Fig. Wc)] and find that 
there is a negative eigenvalue among an infinite number of 
positive ones, which means that there is a single path in the 
configuration space along which the energy decreases. That 
direction belongs to a in = perturbation which takes the 
vortex out of the condensate \^. Nevertheless, as in this case 
there exist no complex eigenvalues, we conclude that the life- 
time of the configuration is only limited by the presence and 
amount of the losses. This is confirmed when take a (1, 0) con- 
figuration, perturb it and study its real time evolution [Fig. 
hid )], and it is indeed consistent with the experiments of Ref. 

We must remark that the existence of a negative eigen- 
value in the spectrum around the (1,0) vortex contradicts 
the belief that the second component could act as a pinning 
potential. Indeed, based on further study we can state that a 
(1,0) vortex remains energetically unstable even for different 
proportions of each specie, from A^i ~ to A^2 — [[13|. 
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FIG. 1. (a) Energy of the ground state, (0,0), (b) energy 
of the (1,0) (solid line) and (0, 1) (dashed line) configurations 
relative to the ground state and (c) frequencies of the linear 
modes around the (1, 0) stable state, with respect to the non- 
linearity, given by g [Eq. (H)] . (d) Upper view of component 
1 1) with a vortex after a strong perturbation of the cloud (see 
above) for g = 1500. Pictures are equispaced by 10 time units 
and each one is 18 x 18 spatial units large. 

Finally we focus on the (0,1) family [Fig. M(a)]. Here 
we find an infinite number of modes with positive energy, 
plus a pair of them with negative frequency and complex 
eigenvalues, Au. Qualitatively, the shape and frequencies of 
the unstable modes are similar to those of the energy de- 
creasing modes of the (1,0) -that is, they are perturbations 



which push the vortex out of both clouds. The difference is 
that due to the imaginary part of those eigenvalues, which is 
Im(A„) = C(0.04) [Fig. |(b)], vortices with unit charge in |2) 
are unstable under a generic perturbation of the initial data. 
This is consistent with the JILA experiments, where the a 
vortex hosted inside a |2) specie was found to be unstable. 

Although, as we mentioned above, the previous result does 
not depend drastically on A'^i being equal to A'^2, we have also 
found that for N2 ^ A'^i the vortex in |2) becomes practically 
stable |l3] -that is, the lifetime is too long to be observed 
in numerical simulations. This is consistent with the limit 
of a single-specie condensate where the unit-charge vortex is 
found to be stable Eol for any scattering length. 
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FIG. 2. (a) Linear response frequencies around the (0, 1) 
unstable state and (b) imaginary part of the unstable mode 
arising from the linear stability analysis, both with respect to 
the nonlinearity g [Eq. (H)] (c-d) Upper view of an unstable 
vortex after a finite perturbation for g — 1500. Shown are 
species |1) (c) and |2) (d). Pictures are equispaced by 10 time 
units and each one is 18 x 18 units large. 

Does the vortex break?.- The linear stability analysis cannot 
be used to draw conclusions about the behavior of the vor- 
tex far from the limit of infinitesimal perturbations. To get 
further insight on the dynamics we have performed a set of nu- 
merical simulations in which we reproduce different realistic 
perturbations of each state. These perturbations vary from 
small displacements of the vortex core which should agree in 
detail with the linear stability results, to strong perturbations 
of the initial data. In particular, we have systematically used 
the procedure of creating a stationary state and then reducing 
the size of the trap, thus transfering the system to a non sta- 
tionary state mM . The numerical simulations have been done 
using a symmetrized split-step Fourier pseudospectral method 
on grids of sizes ranging from 32 x 32 x 32 to 64 x 64 x 64 for 
the three-dimensional setups, and 64 x 64 to 256 x 256 for the 
two-dimensional model. All results were tested on different 
grid sizes and changing time steps to ensure their validity. 

From our numerical simulations we extract several conclu- 
sions. First, the linearly stable state, (1,0), is also a robust 



one and survives to an ample range of perturbations, suffer- 
ing at most a precession of the vortex core plus changes of 
the shapes of both components [Fig. Wd)]. This behavior 
is equally reproduced in two- and three-dimensional simula- 
tions. 

Next we have the unstable states, (0, 1) subject to weak per- 
turbations in either two- or three-dimensional models. In that 
case the unstable configuration develops a simple recurrent 
dynamics, which is well represented by Figs. 0(d-e). There 
we see a first stage where the first component [Fig. 0(d)] 
and the vortex [Fig. 0(e)] oscillate synchronously (the hole 
in 1 2) pins the peak of |1)). These oscillations grow in ampli- 
tude until the linked system spirals out and forms a ying-yang 
which rotates clockwise, one specie chasing the other. Finally 
the first component develops a tail and later a hole which 
traps the second component. That hole is a vortex, which 
somehow has been transfered from |2) to |1). Though it is 
not completely periodic, this mechanism does exhibit some 
recurrence and the vortex returns to |2). 

The preceding behavior persists even for strong perturba- 
tions in a two-dimensional condensate. However, when one 
considers a three dimensional condensate and applies large 
perturbations on the initial data such as the one described 
above, it is possible to find a richer dynamics which include 
the establishment of spatiotemporal chaos. As an example. 
Fig. p^ shows a regime in which more than one vortex are in- 
troduced into the first component for long times. Intuitively, 
this turbulent behavior has two causes. First, due to ener- 
getic considerations there is a bigger overlap of both species 
than in the two dimensional model, as is apparent from the 
pictures. And second and most important, in the a three di- 
mensional environment the first component is more reluctant 
to be dragged by the weak vortex-line from the second com- 
ponent. Thus, as the second component spirals around, it 
is able to shake the first fluid and produce pairs of vortices, 
much like what happens in laser-stirred condensates [h|. 



m% I f ^ <i 



• •» c» M % t 



FIG. 3. Evolution of an unstable (0, 1) configuration sub- 
ject to a finite perturbation after a 50 units run for g — 1500. 
Spatiotemporal chaos with several vortices develops for long 
enough times. Shown are upper views equispaced by 10 time 
units, each one 18 x 18 spatial units large. 

Either phenomena, the transfer of the vortex and the 
chaotic evolution do not not contradict any topological rule or 
conservation law. In fact the angular momentum of each com- 
ponent is no longer a conserved quantity, and the topological 
charge of each specie needs not survive through the evolution. 
Instead, what is conserved is the total angular momentum 
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In Fig. Hwe plot the evolution of the total angular momentum 
and that of each component L^ for a recurrent situation and 
for the case of a spatiotemporally chaotic state [Fig. ra . It can 
be seen that there is a complex interchange of angular mo- 
mentum between both components. The intermediate states 
are topologically nontrivial ones since the phase singularity is 
being transfered from one component to the other. A more 
detailed analysis of this process will be presented elsewhere 
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FIG. 4. Evolution of total angular momenta L^ (solid line, 
well conserved by our numerical scheme), and partial angu- 
lar momenta Li (dotted line) and Li (dashed line), with 
respect to time in adimensional units. Plot (a) belongs to 
a situation with a simple, recurrent transfer of the vortex, 
while (b) corresponds to a chaotic spatiotemporal dynamics 
[Fig. |. 

Conclusions and discussion.- We have analyzed the stabil- 
ity of vortices in multi-component atomic Bose-Einstein con- 
densates, using both two-dimensional and three-dimensional 
sets of coupled Gross-Pitaevskii equations. We prove that 
a vortex in a |1) state is a dynamically stable object even 
though it is not a global minimum of the energy. In contrast, 
we demonstrate that a vortex in the |2) specie is dynamically 
and energetically unstable and tends to spiral out of the con- 
densate. Besides, since our model does not involve dissipative 
effects or asymmetries, one cannot get rid of the vortex angu- 
lar momentum. This leads to a complex dynamics in which 
the vortex is transfered to the |1) state and eventually goes 
back to 1 2). This and other predictions about the dynamical 
behavior can be checked in current experiments. 

We believe that the simple model presented in this paper 
gives a reasonable explanation of the experiments of [hll based 
only on the nonlinear interactions present in mean field the- 
ories. In fact if the instability found in the experiments were 
due to dissipation through a mechanism similar to that pro- 
posed in Ref. [h7| it would affect both the (1,0) and (0,1) 
type of vortices in a similar way, which is not the case. In 
support to our theory we can see from Fig. 3 of Ref. [hll that 
the vortex does not completely escape but some kind of defect 
is formed in the periphery of the |2) component, a behavior 
similar to that found in our real time simulations [Fig 2(c)]. 

Regarding the time scales, our Zmear stability analysis gives 
hfetimes of about 500 milliseconds, which is two or three times 
larger than what is seen in the experiments. Nevertheless this 
is not significant since further study has revealed that the time 



after which the instability affects the system depends dramat- 
ically on the type and the intensity of the perturbation, which 
the linear stability analysis requires to be small. As such it 
is conceivable that the experimental realization of the (0, 1) 
vortex must break sooner: the whole experimental prepara- 
tion takes the system to a state which differs finitely from the 
stationary configuration and are thus more likely to excite the 
unstable mode, even through the preparation process. 

This work has been partially supported by the DGICYT 
under grant PB96-0534. 
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